home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / rng / rand48.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-05-05  |  3.8 KB  |  143 lines
  1. /* rng/rand48.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <math.h>
  22. #include <stdlib.h>
  23. #include <gsl/gsl_rng.h>
  24.  
  25. /* This is the Unix rand48() generator. The generator returns the
  26.    upper 32 bits from each term of the sequence,
  27.  
  28.    x_{n+1} = (a x_n + c) mod m 
  29.  
  30.    using 48-bit unsigned arithmetic, with a = 0x5DEECE66D , c = 0xB
  31.    and m = 2^48. The seed specifies the upper 32 bits of the initial
  32.    value, x_1, with the lower 16 bits set to 0x330E.
  33.  
  34.    The theoretical value of x_{10001} is 244131582646046.
  35.  
  36.    The period of this generator is ? FIXME (probably around 2^48). */
  37.  
  38. static inline void rand48_advance (void *vstate);
  39. static unsigned long int rand48_get (void *vstate);
  40. static double rand48_get_double (void *vstate);
  41. static void rand48_set (void *state, unsigned long int s);
  42.  
  43. static const unsigned short int a0 = 0xE66D ;
  44. static const unsigned short int a1 = 0xDEEC ;
  45. static const unsigned short int a2 = 0x0005 ;
  46.  
  47. static const unsigned short int c0 = 0x000B ;
  48.  
  49. typedef struct
  50.   {
  51.     unsigned short int x0, x1, x2;
  52.   }
  53. rand48_state_t;
  54.  
  55. static inline void
  56. rand48_advance (void *vstate)
  57. {
  58.   rand48_state_t *state = (rand48_state_t *) vstate;
  59.  
  60.   /* work with unsigned long ints throughout to get correct integer
  61.      promotions of any unsigned short ints */
  62.  
  63.   const unsigned long int x0 = (unsigned long int) state->x0 ;
  64.   const unsigned long int x1 = (unsigned long int) state->x1 ;
  65.   const unsigned long int x2 = (unsigned long int) state->x2 ;
  66.  
  67.   unsigned long int a ;
  68.   
  69.   a = a0 * x0 + c0 ;
  70.   state->x0 = (a & 0xFFFF) ;
  71.  
  72.   a >>= 16 ;
  73.  
  74.   /* although the next line may overflow we only need the top 16 bits
  75.      in the following stage, so it does not matter */
  76.  
  77.   a += a0 * x1 + a1 * x0 ; 
  78.   state->x1 = (a & 0xFFFF) ;
  79.  
  80.   a >>= 16 ;
  81.   a += a0 * x2 + a1 * x1 + a2 * x0 ;
  82.   state->x2 = (a & 0xFFFF) ;
  83. }
  84.  
  85. static unsigned long int 
  86. rand48_get (void *vstate)
  87. {
  88.   unsigned long int x1, x2;
  89.  
  90.   rand48_state_t *state = (rand48_state_t *) vstate;
  91.   rand48_advance (state) ;
  92.  
  93.   x2 = (unsigned long int) state->x2;
  94.   x1 = (unsigned long int) state->x1;
  95.  
  96.   return (x2 << 16) + x1;
  97. }
  98.  
  99. static double
  100. rand48_get_double (void * vstate)
  101. {
  102.   rand48_state_t *state = (rand48_state_t *) vstate;
  103.  
  104.   rand48_advance (state) ;  
  105.  
  106.   return (ldexp((double) state->x2, -16)
  107.       + ldexp((double) state->x1, -32) 
  108.       + ldexp((double) state->x0, -48)) ;
  109. }
  110.  
  111. static void
  112. rand48_set (void *vstate, unsigned long int s)
  113. {
  114.   rand48_state_t *state = (rand48_state_t *) vstate;
  115.  
  116.   if (s == 0)  /* default seed */
  117.     {
  118.       state->x0 = 0x330E ;
  119.       state->x1 = 0xABCD ;
  120.       state->x2 = 0x1234 ;
  121.     }
  122.   else 
  123.     {
  124.       state->x0 = 0x330E ;
  125.       state->x1 = s & 0xFFFF ;
  126.       state->x2 = (s >> 16) & 0xFFFF ;
  127.     }
  128.  
  129.   return;
  130. }
  131.  
  132. static const gsl_rng_type rand48_type =
  133. {"rand48",            /* name */
  134.  0xffffffffUL,            /* RAND_MAX */
  135.  0,                /* RAND_MIN */
  136.  sizeof (rand48_state_t),
  137.  &rand48_set,
  138.  &rand48_get,
  139.  &rand48_get_double
  140. };
  141.  
  142. const gsl_rng_type *gsl_rng_rand48 = &rand48_type;
  143.